
function dxx = Cr3bp_Stm2D(tau , xx , mu)
%
% CR3BP动力学模型，状态+STM，2D
%
% 例子：
% aux = crtbpAuxEM;
% crtbpMarkEM
% xx0 = [0.75 , 0 , 0.2 , 0.15 , 0 , 0]';
% yy0 = [xx0 ; reshape(eye(6) , 36 , 1)];
% options = odeset('Reltol' , 1e-13, 'AbsTol' , 1e-13);
% [tt , yy] = ode113(@crtbpEqmSTM3D , [0 , 1e2] , yy0 , options , aux.mu);
% plot3D(yy(: , 1 : 3) , 'b' , 1);
%
% 时间：2024年06月03日
% 作者：张晨
% 单位：CSU
% 邮箱：chenzhang@csu.ac.cn
%%%%%%%%%%%%%%%%%%%%%%%

% % 积分时间太长，跳过
% if nargin == 4
%     ode_t = clock;
%     if etime(ode_t , ode_t0) > 1
%         dxx = [];
%         return
%     end
% end

% % 配置参数
% mu = aux.ems.mu;
% LU = aux.ems.LU;
% P1_req = aux.planet.req(3);
% P2_req = aux.planet.req(10);

% 提取变量
x = xx(1);
y = xx(2);
vx = xx(3);
vy = xx(4);

% % P1/P2碰撞检测
% if aux.impact_IO == 1
%     impactP1 = sqrt( (x + mu)^2 + y^2 ) - P1_req / LU;
%     impactP2 = sqrt( (x + mu - 1)^2 + y^2 ) - P2_req / LU;
%     if impactP1 * impactP2 < 0
%         dxx = [];
%         return;
%     end
% end

% 动力学方程
r1_3pow = ((x + mu)^2 + y^2)^(1.5);
r2_3pow = ((x + mu - 1)^2 + y^2)^(1.5);
r1_5pow = r1_3pow^(5/3);
r2_5pow = r2_3pow^(5/3);

G11 = 1 - (1 - mu) / r1_3pow + 3 * (1 - mu) * (x + mu)^2 / r1_5pow - mu / r2_3pow + 3 * mu * (x + mu - 1)^2 / r2_5pow;
G12 = 3 * (1 - mu) * (x + mu) * y / r1_5pow + 3 * mu * (x + mu - 1) * y / r2_5pow;
G22 = 1 - (1 - mu) / r1_3pow + 3 * (1 - mu) * y^2 / r1_5pow - mu / r2_3pow + 3 * mu * y^2 / r2_5pow;

Df = [0 , 0 , 1 , 0;
    0 , 0 , 0 , 1;
    G11 , G12 , 0 , 2;
    G12 , G22 , -2 , 0];
Phi = reshape(xx(5 : 20) , 4 , 4);
dPhi = Df * Phi;

dxx = zeros(20 , 1);
dxx(1 : 4) = [vx;
    vy;
    x - (1 - mu) * (x + mu) / r1_3pow - mu * (x + mu - 1) / r2_3pow + 2 * vy;
    y - (1 - mu) * y / r1_3pow - mu * y / r2_3pow - 2 * vx];
dxx(5 : 20) = dPhi(:);

end
